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One of the most important problems in magnetically confined fusion research is the turbulent 
transport of heat observed in present-day Tokamak experiments. Plasma turbulence is very 
challenging from a theoretical point of view due to the nonlinearity and high dimensionality of the 
governing equations. Likewise, experimental measurements in the core region of a Tokamak are 
limited by the extremely high temperatures -10^ °K. The level of both theoretical and empirical 
difficulty highlights the important role numerical simulations can potentially play in developing a 
predictive model for turbulent heat transport. Such a model would dramatically reduce 
uncertainties in design and may lead to enhanced operating regimes, both of which would reduce the 
i and expense of a fusion reactor. Simulating this highly non-linear behavior requires solving for 

^*H?ie perturbations of the phase space distribution function in five dimensions (via Sf methods using 
gyrokinetic formalism). We use a particle-in-cell approach to solve the equations. The code has 
^Jt^een parallelized for a variety of architectures using a 1-D domain decomposition along the toroidal 
^axis, for which the number of particles on each processor remains approximately constant — 
1 minimizing load imbalances and interprocessor communication — making this problem ideally suited 
kJo massively parallel architectures. We present the performance of the program for different 
^numbers of processors and problem sizes. 

O 

^jeometry 

^ The essential geometry of a Tokamak is that of a torus defined by a major radius R and a minor 
Ci»dius a (see Figure 1). The ions within the plasma move rapidly around the torus, gyrating 
tightly along the magnetic field lines, like rings on a wire. The radius of gyration p is set by the 
velocity, mass and charge of the particle. The essential scale of the system is set by a/p. Typically 
l/^ ~ 270 cm, a - 85 cm, and p ~ 0.15 cm. 

^s. \ The simplest arrangement of field lines is obtained by wrapping current carrying wires tightly 
\ jlround the minor axis of the torus creating straight magnetic field lines aligned with the torus. 
N^filnfortunately, the magnetic field exerts a greater force on the inside of the torus, which causes the 


0ns in the plasma to drift across the field lines. This problem can be partially alleviated by twisting 
^Jjie field lines into a helical shape so that the drifts approximately cancel. A byproduct of this 
^twisting is that the trajectories of particles are now far more complex and are susceptible to a wide 
V^ange of instabilities, which tend to grow along the toroidal modes of the Tokamak. Figure 2 
C^Shows a simplified example of a typical toroidal mode. This could be, for example, the linear 
"-■o^rowth phase of an ion density perturbation. Generally, the isosurfaces of a perturbation follow the 
Cyftelical shape of the magnetic field lines. 

o 


implementation 


Particle-in-cell codes have been used in the plasma physics community for several decades. The 
^^Sssence of the computational problem involves solving for the trajectories of particles with mass m 
f anil charge q using equations analogous to 


dVi(t)/dt = (q/m) E(Xj(t)) , dX](t)/dt = Vi(t) , 

V-E = 4jc q n(x) , n(x) = Zi S(x - Xi) , 

where the subscript i refers to the ith particle and the electric field at the particle's position is 
E(xj(t)) is found by interpolating from fields previously calculated on a grid. The interpolation is a 
gather operation which involves indirect addressing and a substantial part of the computation time is 
spent here. The fields created by the particles are found from Poisson's equation. This equation is 
solved on a grid, using Fourier transform methods. Typically, the time for the field solver is not 
large. The source term n(x) is calculated from the particle position by an inverse interpolation 
where S is a particle shape function. This is a scatter operation which also involves indirect 
addressing and consumes a substantial part of the computation time. 



Figure 1. Schematic drawing of a Tokamak with major radius R 
and minor radius a, and the gyrating path of a particle (black line) 
along a magnetic field line (gray line) with a gyration radius p. 



Figure 2. Example of a m=3, n=2 mode in a torus, m is the 
poloidal mode number and n is the toroidal mode number. The 
isosurface shows the geometry of a typical growing standing wave 
within the plasma, it also traces the magnetic field lines. 



Figure 3. Ion density in the linear phase of the simulation showing 
the ion-temperature-gradient instability. The mode is very elongated 
in the direction following the magnetic field lines. 


Results 


When the simulations start with a very small initial perturbation, there are generally two 
identifiable phases of the run. First is the linear phase, where modes (i.e., standing waves) grow 
exponentially. Linear modes can also be found using lower-dimensional, time-independent 
eigenvalue techniques and are fairly well understood theoretically. The second phase is the turbulent 
stationary state, during which the growth of the dominant linear modes saturates and the system 
settles down to a statistical steady state. The ion density in the linear mode is shown on surfaces 
throughout the Tokamak in Figure 3. 

In our implementation, a ID domain decomposition along the toroidal axis is used, which is 
significantly easier to program than a full 3D decomposition. In addition, along this axis the 
number of particles per processor remains relatively constant, which minimizes load imbalance. 
The most recent efforts have been in porting the code to the Cray T3D using FortranW and the 
PVM message passing library. Production runs on the T3D show a performance of approximately 
14.4 MFlops per processor (10% of the theoretical peak of 150 MFlops, which is typical of these 
types applications). To test the scalability of the code, 12 runs were timed with different numbers 
of processors and problem sizes. To first order, the problem size is given by the total number of 
particles. The total amount of computer resource consumed is the time per step multiplied by the 
number of processors. In a perfectly scalable code, the resource consumption should be 
proportional to the problem size, which is shown by the straight line in Figure 4. 


Figure 4. Total processor time used vs. problem size. The 
straight line indicates the expected time from scaling up the smallest 
simulation. 


